title: “Week 3: Use R as GIS”

output: github_document

knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/',  warning=FALSE, message=FALSE)

## [[1]]
## [1] "sp"        "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[2]]
## [1] "rgdal"     "sp"        "stats"     "graphics"  "grDevices" "utils"    
## [7] "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "rgeos"     "rgdal"     "sp"        "stats"     "graphics" 
##  [6] "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[4]]
##  [1] "maptools"  "rgeos"     "rgdal"     "sp"        "stats"    
##  [6] "graphics"  "grDevices" "utils"     "datasets"  "methods"  
## [11] "base"     
## 
## [[5]]
##  [1] "classInt"  "maptools"  "rgeos"     "rgdal"     "sp"       
##  [6] "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [11] "methods"   "base"     
## 
## [[6]]
##  [1] "RColorBrewer" "classInt"     "maptools"     "rgeos"       
##  [5] "rgdal"        "sp"           "stats"        "graphics"    
##  [9] "grDevices"    "utils"        "datasets"     "methods"     
## [13] "base"        
## 
## [[7]]
##  [1] "GISTools"     "MASS"         "RColorBrewer" "classInt"    
##  [5] "maptools"     "rgeos"        "rgdal"        "sp"          
##  [9] "stats"        "graphics"     "grDevices"    "utils"       
## [13] "datasets"     "methods"      "base"        
## 
## [[8]]
##  [1] "maps"         "GISTools"     "MASS"         "RColorBrewer"
##  [5] "classInt"     "maptools"     "rgeos"        "rgdal"       
##  [9] "sp"           "stats"        "graphics"     "grDevices"   
## [13] "utils"        "datasets"     "methods"      "base"        
## 
## [[9]]
##  [1] "raster"       "maps"         "GISTools"     "MASS"        
##  [5] "RColorBrewer" "classInt"     "maptools"     "rgeos"       
##  [9] "rgdal"        "sp"           "stats"        "graphics"    
## [13] "grDevices"    "utils"        "datasets"     "methods"     
## [17] "base"        
## 
## [[10]]
##  [1] "ggmap"        "ggplot2"      "raster"       "maps"        
##  [5] "GISTools"     "MASS"         "RColorBrewer" "classInt"    
##  [9] "maptools"     "rgeos"        "rgdal"        "sp"          
## [13] "stats"        "graphics"     "grDevices"    "utils"       
## [17] "datasets"     "methods"      "base"

Spatial Objects

	Without attributes	With attributes
Points	SpatialPoints	SpatialPointsDataFrame
Lines	SpatialLines	SpatialLinesDataFrame
Polygons	SpatialPolygons	SpatialPolygonsDataFrame
Raster	SpatialGrid	SpatialGridDataFrame
Raster	SpatialPixels	SpatialPixelsDataFrame

LubbockBlock<-readShapePoly("Data/LubbockBlockNew.shp") #read polygon shapefile
class(LubbockBlock)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

HouseLocation<-read.csv("Data/HouseLatLon.csv") #read GPS data
class(HouseLocation)

## [1] "data.frame"

coordinates(HouseLocation)<-c('Lon', 'Lat')
class(HouseLocation)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

cropland<-raster("Data/Lubbock_CDL_2013_USDA.tif")
class(cropland)

## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"

tmin <- getData("worldclim", var = "tmin", res = 10)  # this will download 
class(tmin)

## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

LubbockBlock<-readOGR("./Data", "LubbockBlockNew") #read polygon shapefile

## OGR data source with driver: ESRI Shapefile 
## Source: "./Data", layer: "LubbockBlockNew"
## with 167 features
## It has 69 fields

class(LubbockBlock)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

Mapping with R

Basic Mapping

LubbockBlock<-readShapePoly("Data/LubbockBlockNew.shp") #read polygon shapefile
plot(LubbockBlock,axes=TRUE, col=alpha("gray70", 0.6)) #plot Lubbock block shapefile
#add title, scalebar, north arrow, and legend
HouseLocation<-read.csv("Data/HouseLatLon.csv") #read GPS data
price<-HouseLocation$TotalPrice
nclr<-5
priceclr<-brewer.pal(nclr, "Reds")
class<-classIntervals(price, nclr, style="quantile")
clocode<-findColours(class, priceclr)

points(HouseLocation$Lon, HouseLocation$Lat, pch=19, col=clocode, cex=0.5) #add houses on top of Lubbock block shapefile
title(main="Houses on Sale in Lubbock, 2014") 

legend(-101.95, 33.65, legend=names(attr(clocode, "table")), fill =attr(clocode, "palette"), cex=0.5, bty="n")
#map.scale(x=-101.85, y=33.49,0.001,"Miles",4,0.5,sfcol='red')
north.arrow(xb=-101.95, yb=33.65, len=0.005, lab="N")

#plot raster
plot(cropland)

#plot raster stack
tmin <- getData("worldclim", var = "tmin", res = 10)  # this will download 
plot(tmin)

Mapping with static Google Maps

Mapping with dynamic Google Maps

Changing map projections

#project a vector 

boudary=readShapePoly('Data/boundary');
proj4string(boudary) <-CRS("+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
proj4string(boudary)

## [1] "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

boudaryProj<-spTransform(boudary, CRS("+init=epsg:3857"))
proj4string(boudaryProj)

## [1] "+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"

#project a raster
proj4string(cropland)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

plot(cropland)

aea <- CRS("+init=ESRI:102003")   #Albert equal area
projCropland=projectRaster(cropland, crs=aea)
plot(projCropland)

Spatial analysis with R

vector analysis (overlay)

#project a vector 

# Datasets
#  * CSV table of (fictionalized) brown bear sightings in Alaska, each
#    containing an arbitrary ID and spatial location specified as a
#    lat-lon coordinate pair. 
#  * Polygon shapefile containing the boundaries of US National Parks
#    greater than 100,000 acres in size. 

bears <- read.csv("Data/bear-sightings.csv")
coordinates(bears) <- c("longitude", "latitude")

# read in National Parks polygons
parks <- readOGR("Data", "10m_us_parks_area")

## OGR data source with driver: ESRI Shapefile 
## Source: "Data", layer: "10m_us_parks_area"
## with 61 features
## It has 8 fields

# tell R that bear coordinates are in the same lat/lon reference system as the parks data 
proj4string(bears) <- proj4string(parks)

# combine is.na() with over() to do the containment test; note that we
# need to "demote" parks to a SpatialPolygons object first
inside.park <- !is.na(over(bears, as(parks, "SpatialPolygons")))

# calculate what fraction of sightings were inside a park
mean(inside.park)

## [1] 0.1720648

## [1] 0.1720648

# determine which park contains each sighting and store the park name as an attribute of the bears data
bears$park <- over(bears, parks)$Unit_Name

# draw a map big enough to encompass all points, then add in park boundaries superimposed upon a map of the United States
plot(bears)
map("world", region="usa", add=TRUE)
plot(parks, border="green", add=TRUE)
legend("topright", cex=0.85, c("Bear in park", "Bear not in park", "Park boundary"), pch=c(16, 1, NA), lty=c(NA, NA, 1), col=c("red", "grey", "green"), bty="n")
title(expression(paste(italic("Ursus arctos"), " sightings with respect to national parks"))) 

# plot bear points with separate colors inside and outside of parks
points(bears[!inside.park, ], pch=1, col="gray")
points(bears[inside.park, ], pch=16, col="red")

# write the augmented bears dataset to CSV
write.csv(bears, "bears-by-park.csv", row.names=FALSE)

# ...or create a shapefile from the points
writeOGR(bears, ".", "bears-by-park", driver="ESRI Shapefile")

Raster analysis

tmin=getData('worldclim', var='tmin', res=10)

# Raster calculator
diff=tmin$tmin1 - tmin$tmin10

## the following code is faster for large datasets. 
overlay(tmin$tmin1, tmin$tmin10, fun=function(x,y){return (x-y)})

## class       : RasterLayer 
## dimensions  : 900, 2160, 1944000  (nrow, ncol, ncell)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : -343, 102  (min, max)

elevation <- getData("alt", country = "ESP")
slope <- terrain(elevation, opt = "slope")
aspect <- terrain(elevation, opt = "aspect")
hill <- hillShade(slope, aspect, 40, 270)
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Spain")
plot(elevation, col = rainbow(25, alpha = 0.35), add = TRUE)

#contours

contour(elevation)

#crop raster
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Spain")
plot(elevation, col = rainbow(25, alpha = 0.35), add = TRUE)
extent=drawExtent()
cropElev <- crop(elevation, extent)
plot(cropElev)